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Abstract 

We consider the exponential matrix representing the dynamics of the 
Fermi-Bose model in an undepleted bosonic held approximation. A recent 
application of this model is molecular dimers dissociating into its atomic 
compounds. The problem is solved in D spatial dimensions by dividing 
the system matrix into blocks with generalizations of Hankel matrices, 
here refered to as D-block-Hankel matrices. The method is practically 
useful for treating large systems, i.e. dense computational grids or higher 
spatial dimensions, either on a single standard computer or a cluster. In 
particular the results can be used for studies of three-dimensional phys- 
ical systems of arbitrary geometry. We illustrate the generality of our 
approach by giving numerical results for the dynamics of Glauber type 
atomic pair correlation functions for a non-isotropic three-dimensional 
harmonically trapped molecular Bose-Einstein condensate. 



1 Introduction 

The Fermi-Bose model under study here forms the underlying basis for a range 
of phenomena in condensed matter and ultra-cold atomic physics. It was pro- 
posed in the context of high-temperature superconductivity by Friedberg and 
Lee PQ, and in ultracold gases it corresponds to the theory of resonance su- 
perfluidity with Feshbach molecules (2] El S] - The latter forms the basis of a 
model for describing the physics of the BCS-BEC crossover [5]. More recently, 
the fermion-boson model has been used for analyzing the decay of double oc- 
cupancies (doublons) [5] in a driven Fermi- Hubbard system [TJ. The particu- 
lar situation that we concentrate on in this article corresponds to spontaneous 
dissociation of a Bose-Einstein condensate of molecular dimers into fermionic 
atoms [H [9|. This process represents a fermionic counterpart of parametric 
down-conversion in quantum optics. After recent experimental achivements of 
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molecular dissociation [101 [TO [TO \T3\ , our aim here is to take the theory for 
numerical modeling of dissociation into fermionic atoms from the state of prin- 
cipally possible to the state of being useful in practice. 



1.1 Effective quantum field theory 

Let us here in brief present the three-wave interaction type Hamiltonian of 
interest in this work [TO [TO HU . 

H = H - ih X / (*JM 2 - HH^o) ■ (1) 

Here all the quadratic terms are collected in Hq and contains kinetic- and poten- 
tial energy terms, ^o( x ;*) stands for a bosonic field operator, whereas ^(x, t) 
(j = 1,2) describe two particle fields that can be two fermions (bosons) in dif- 
ferent spin states, finally \ is the strength of the fermion-boson (boson-boson) 
coupling term. 



1.2 Applications of the Fermi-Bose model 

In modern condensed matter physics the Fermi-Bose model have two major areas 
of applicability. First the so called "s-channel" model in high-temperature su- 
perconductivity [T] . In this context it model the formation dynamics of bosonic 
Cooper-pairs, 

i — H — * 

Cooper-pair 

where the two atomic particles, the electrons, are fermions. 
In the field of ultra-cold atomic physics, it can model the dissociation of ultra- 
cold bosonic molecules [TO [18] , 




ATOM 1 



MOLECULE 



ATOM 2 



and hence we allow here for the atomic particles to be either two fermions 
or two bosons. 



1.3 Computational methods for large systems 

Various generalizations of time-dependent DMRG PH] to higher dimensional 
systems is a topic of large present interest in the computational physics commu- 
nity [20], but have not yet reached a useful status for large higher dimensional 
systems as of interest here. 
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Methods for bosonic evolution based on stochastic differential equations 
(SDEs), with the ability for independent stochastic trajectories to be carried 
out on computer clusters for large systems, are succesful in many situations 
[TT1 [T51 [22 [221 [231 [2H 121] but are restricted to simulations of molecules dissoci- 
ating into bosonic atoms. 

In this article we outline a method to study the fermionic time evolu- 
tion for effective Heisenberg equations that are linear in creation and annihila- 
tion operators. We apply the method to the evaluation of analytic short-time 
asymptots for the Glauber's second order correlation functions [23], for a non- 
isotropic three-dimensional molecular condensate dissociating into fermionic 
atoms, against numerical data. We focus first on a general formulation in a 
D-dimensional Cartesian momentum base, such that convenient numerical so- 
lutions can be directly obtained for an arbitrary shaped bosonic field. However, 
for a specific application the numerical performance may be further increased 
by using additional geometrical symmetries or by formulating the equations 
in a different basis. Except for generating valuable results in certain physical 
regimes, the method presented here can also be useful as a reference for vali- 
dation of approximate analytic results or to evaluate more advanced numerical 
methods, such as for example the Gaussian fermionic phase-space representa- 
tion (GPSR) that have recently been applied to the Fermi-Bose model with a 
uniform bosonic field [251 [23 and an implementation of GPSR for dissociation 
from a non-uniform molecular BEC is in progress. Turning off stochastic terms 
in the SDEs of the GPSR, one obtain in effect the so called pairing mean-field 
theory (PMFT) [H [26l [28]. Furthermore, keeping the molecular variables in 
PMFT undepleted will be equivalent to the formulation in the present article 
but is numerically less suitable. 

There is an obvious computational advantage to solve for the single operator 
dynamics, instead of solving for pairs of operators, the latter is done e.g. in the 
PMFT (and GPSR) discussed above. To give one relevant example, to simulate 
a three-dimensional non-uniform field on a Cartesian momentum grid of size 
100 x 100 x 100 requires in effect, as we will show in this paper, only to be able to 
store part of a (sparse) n x n D-block-Hankel matrix of the size n = 100 3 = 10 6 
and to multiply this D-block-Hankel matrix with vectors only for each final 
time-point of interest. On the other hand c-number based mean-field methods 
for Fermi-Bose systems like PMFT, where the basic variables represents pairs of 
operators with two indices, requires in this case to propagate ~ 10 12 variables 
in time through sufficiently many small time-steps up to the final time-point 
of interest. The methods for stochastic evolution that represents single bosonic 
field operators with complex stochastic fields mentioned above, see [21] for a 
recent review, do not have any direct corresponding useful method for fermions. 
For example the GPSR that can treat the quantum dynamics of the Fermi-Bose 
model exact [2"r?l [771 HZ5] , involves basis elements that represents pairs of single 
operators [23 123 12S1 [33 1311 132] > hence it is generally more restrictive in the 
size of the computational grid. 

Due to the existence of several related methods for bosonic dynamics, we fo- 
cus in the present article in particular on applications to systems with fermionic 



3 



atomic operators. However, the specific analogue theory of two distinguishable 
bosonic atomic operators is also presented for comparison, since its formulation 
differ only with a sign. 

Advantages with the method presented here are i) to be able to study effects 
of quantum statistics [by chosing q = — 1 (q = 1) for fermionic (bosonic) atomic 
particles]; ii) to obtain accurate correlation functions for short times i.e. for a 
small number of atomic particles, where stochastic evolution methods have a low 
signal to noise ratio; and Hi) to conviniently treat systems with a moderate size 
(relative to the RAM memory) of the computational grid (e.g. lower dimensional 
systems), where deterministic values of the observables are obtained fast even 
on a single standard PC. 

The article is organized as follows. Section [5] provides the Heisenberg equa- 
tion of motion to be treated within the concept of molecular dissociation into 
fermionic (bosonic) atoms and also briefly mention dissociation into two indis- 
tinguishable bosonic atoms and the related problem of condensate collision. In 
section |3] we relate blocks in the system matrix responsible for atom-molecule 
coupling to the so called £>-block-Hankel matrices. In section 2] we formu- 
late and prove results for the structure of the solution of Heisenberg equation 
of motion. Section [5] provides the reader with the practical details of how 
to use the mathematical results in obtaining physical observables. We apply 
the general method to a three-dimensional non-isotropic harmonically trapped 
Bose- Einstein condensate in section [B] and present numerical results for atomic 
correlation functions that we compare to recently derived analytic short-time 
asymptotes. Finally the article is summarised in section 

2 The momentum-space operator equations 

To model the dissociations of a Bose-Einstein condensate of diatomic molecules 
into pairs of constituent atoms, we start with the following effective quantum 
field theory Hamiltonian 

J ^=0,1,2 lm i 

-ift X (*S*i*2-*£*l*o)}- (2) 

Here we assume that the molecules are made of either two distinguishable 
bosonic atoms or two fermionic atoms in different spin states. In both cases, 
\E'o( x , t) is a bosonic field operator for the molecules, satisfying the standard 
commutation relations [^o( x , t), ^ (x.',t)} = 5 D (x — x'), with D being the spa- 
tial dimension of the system (33]. The atomic field operators, ^(x, t) (j = 1,2), 
satisfy either bosonic commutation or fermionic anti-commutation relations, 
[§;(x,i),§}(x',i)] = %<5 D (x-x') and [$J(x,t),$J(x',t)] = [*i(x,t), §,(x',i)] = 

Oor{* i (x,0,*}(x',t)} = a«tf fl (x-x')a^ 
0, depending on the underlying statistics |33) . 
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The first term in the Hamiltonian © describes the kinetic energy where the 
atomic masses are toi and TO2, whereas the molecular mass is too = mi + 771.2. 
For simplicity, we will consider the case of equal atomic masses (same isotope 
atoms), with mi = to 2 = m a and too = 2m a . 

The coupling constant x = Xd is responsible for coherent conversion of 
molecules into atom pairs, e.g. via optical Raman transitions, an rf transition, 
or a Feshbach resonance sweep and microscopic expressions for x can be found in 
|28| and references therein. The detuning f2 is defined to give the overall energy 
mismatch 2ft£l between the free two-atom state in the dissociation threshold and 
the bound molecular state (including the relative frequencies of the Raman lasers 
or the frequency of the rf field, again see [35J and references therein for details). 
Unstable molecules, spontaneously dissociating into pairs of constituent atoms, 
correspond to ft < 0, with 2ft|Sl| being the total dissociation energy. 

In what follows we will treat the dissociation dynamics in the undepleted 
molecular condensate approximation in which the molecules are represented as 
a fixed classical field. The approximation is valid for short enough dissociation 
times during which the converted fraction of molecules does not exceed about 
10% (24j [28j [34]. In this regime the dissociation typically produces low den- 
sity atomic clouds for which the atom-atom s-wave scattering interactions are 
negligible. For dissociation into bosonic atoms, also the effects of atom-atom 
s-wave scattering have been investigated in |23l |2~41 |2"5] , and the validity of the 
undepleted molecular condensate approximation was found to still hold for a 
converted fraction of molecules up to about 5%. Even though 5%— 10% conver- 
sion efficiencies seem small, nevertheless they can produce mesoscopic ensem- 
bles of pair-correlated atoms with interesting quantum statistics and nontrivial 
many-body correlations if one starts with large-enough molecular condensates, 
such as containing at least 10 4 — 10 5 molecules. Further on, for dissociation 
into fermionic atoms in the regime of Pauli-blocking, where the atomic occu- 
pation numbers are strictly limited by the number of available states [28], the 
undepleted field approximation is expected to be accurate even for large times 
[HI HH1 HE] • In addition to the issue about temporal depletion, the atom- molecule 
interactions will initially appear as an effective spatially dependent detuning due 
to the mean- field interaction energy [35] . To neglect this effect can be motivated 
by operating at relatively large absolute values of the dissociation detuning |f2| 
so that it dominates the mean-field energy shift |24[ 134]. 

The trapping potential for preparing the initial molecular BEC - with any 
residual atoms being removed - is omitted from the Hamiltonian ([2]) since we 
assume that once the dissociation is invoked, the trapping potential is switched 
off, so that the dynamics of dissociation is taking place in free D-dimensional 
space. We assume that the switching on of the atom-molecule coupling x an d 
switching off of the trapping potential is done in a sudden jump at time zero. 
Accordingly the preparation stage is reduced to assuming a certain initial state 
of the molecular BEC in a trap, after which the dynamics is governed by the 
Hamiltonian @. 



2.1 Heisenberg equations in the undepleted molecular con- 
densate approximation 

From the Heisenberg equation of motion with the Hamiltonian taken from ([2]) . 
we have for the three field operators 



(x) 

at 



*i (x) , H 



3 = 0, 1, 2. 



(3) 



In order to obtain linear operator equations the undepleted molecular field 
approximation is first invoked as follows. Assuming that the molecules are 
in a coherent state initially, the density profile po ( x ) is in principal given by 
the ground state solution of the standard Gross-Pitaevskii equation. We then 
replace the molecular field operator by its coherent mean-field complex function 
[23] , the so called condensate wave-function, 



o(x, t) -> (tf (x, t)> = * (x, 0) = vWWexp (i0 (x)) . 



(4) 



From @, and Q we then write down the Heisenberg equations for the 
remaining two coupled atomic field operators as follows 



a*i(x,t) 
at 



2m a 



n 



$i (x, t) + qxVPo (x) exp (%0 (x)) & 2 (x, t) , (5) 



dt 



2m a 



n 



% (x, t) + X VMx) cxp (x)) (x, t) . (6) 



The sign given by q in the second term in (JS|) is q = — 1 for fermionic and q = 1 
for bosonic atoms throughout the paper, as a consequence of different operator 
(anti-) commutator relations. Multiplying |[5} and © with L~ D / 2 exp(— ik ■ x), 
where V = L D is the quantization volume and we assume for simplicity L to be 
the spatial length of the system in any direction, followed by integration over x, 
we can interpret the differential equations for a given k in terms of the Fourier 
operators in momentum space 



1 



L D / 2 



dx ^ (x, t) exp(— ik • x). 



(7) 



The operators Ok,j(t) satisfy the usual (commutation-) anti-commutation rela- 



= <5?j<$k.k' and 



t k,i' "k'.jJ 



[Ok.ij °k' 



(i.e. for 



tions [ok,( 

q = —I, [ , ] +1 = {,}). Since the effective Hamiltonian corresponding to ([5])- 
([6]) is quadratic in the field operators, higher-order moments or expectation 
values of products of creation and annihilation operators will factorize accord- 
ing to Wick's theorem into products of the normal and anomalous densities 

"■k.k'j ,j) an d m k,k' = (&k,i^k',2)- Upon applying the Fourier trans- 

form, the equations ([5]) and (j6]) become 
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daw. 



- = -iAkflk,i + (?k y^Sk'+kO^, 2 , (8) 



at ^-^ 

k 



da} 

- K 5k'+kflkM + iAka^. (9) 



dt 



k' 



where * denotes the complex conjugate. The kinetic part is = fl+h |k| 2 / (2m Q ). 
The effective atom-molecule coupling constant is k = X/L°/ 2 EH]. Finally the 
complex Fourier coefficients of the condensate wave function describing the 
molecular mean-field in momentum-space is defined analogue to ([7]) 

■9 k = JW/2 J dx ^ P° ( x ) exp ( iS ( x ) ~ ik ■ x ) ' ( 10 ) 

For a real condensate wavefunction we have 5^'+k = ff— (k'+k)j while a non-zero 
phase-function 9 (x) can be used to represent for example an initial vortex state 
as in [36]. We give the general theory for a static complex condensate wave- 
function = v Po ( x ) e x P ( x )) m the following, while we for the numerical 
example in section [B] choose = \J Po ( x ) to be real. For the case where the 
molecular density is uniform and constant po: there is only one non-zero Fourier 
coefficient and only operators with index k and k' = k couples in ([H])-©. 
Finally, if the molecular density is non-uniform and time dependent \1/ (x, t) = 
y/p(x.,t) exp (i9 (x, f ) ) , we need to determine the Fourier coefficients ([TO)) for 
each time-step in an iterative process where the dynamics of the molecular 
mean-field is taken into account. This is a topic for future work, however, the 
major effects to be taken into account are: expansion of the condensate |37) ; 
temporal depletion of the number of molecules [23l [24j [25] [26] [28] ; and finally 
the spatial dependences, such as a larger local dissociation rate in the points 
of highest initial molecular density [25]. Defining such an iterative process, the 
results of the present article will be applicable in each time-step of the evolution. 

When the molecular condensate consists of pairs of indistinguishable bosonic 
atoms of a single spin-state, for example 87 Rb2 [11], Heisenberg equation for 
one bosonic operator describe the dissociation dynamics [231 [28] . A similar 
set of Heisenberg equations for bosonic operators of one spin-state can also be 
formulated for the dynamics of condensate collisions within the time-dependent 
Bogoliubov approach [211 12"21 13"51 139| . Hence, also for this problem the method 
discussed in the present article can be applied |40| . 



2.2 Uniform molecular field 

For a reference we start with presenting the relevant results for a size-matched 
uniform system. Size-matched here means that the spatial dimensions of the 
uniform bosonic field L eq , with the same (central-) particle density po, are chosen 
such that the initial number of molecules are the same as for the non-uniform 
system of interest, hence N mo i = poL® q . 
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Under the condition of a uniform molecular field, i.e. that do not depend 
on the spatial coordinates "Jo (x) = \&o = ^/f>o, equations (J5J)- ([9]) with ini- 
tial vacuum states for the atoms have analytical solutions for the normal- and 
anomalous atomic moments n kcr = = ^a kl flk.i^ = (tf-u i^-k,i 

(a k2 ak,2y = (ol kj2 o_k,2) respectively m k = (ok,iO_k,2) = (a-k,i3k,2)- The 
related (PMFT) complex differential equations with initial conditions nk, CT (0) = 
m k (0) = are [28] 



dt 

dmir 



2 5o Re{m k }, (11) 



2iA k m k + g Q (1 + q2n Ka ) , (12) 

where go = ngo = Xy/Po- The corresponding solutions to (jll|) - (fT2|) . calculated 
explicitly in section f5. 21 are 

9o -2 



"k,. = A2 nn 2 sin ( V A k - Ml t ) , (13) 
^k 19o 



mk = y/sJ°- qgl c ° s yJ A *~ q9 ° 7 sin (y/ A l-wo 1 

~Mf^ sin2 G /x ^ t )- (14) 

Note that for bosons (q = 1), e.g.. the resonance mode (A k = 0) leads to a 
Bose-enhancement effect in the atomic occupations, described by n ko (T (t) = 
sinh 2 (got) which grows exponentially with time, consequently this illustrate 
that the undepleted molecular field approximation (Q| is only realistic for a short 
time. In contrast to this, for fermionic atoms, the atomic occupations undergo 
sinusoidal oscillations and can be kept to a small fraction of the number of 
molecules also for large times. As noted in [9] the moments (fT3|) and (fT4|) fulfills 
the equality 



|m k | 2 = (1 + qn^.v) , (15) 



for a uniform molecular field. 



3 Non-uniform molecular field 

We here give the theoretical framework needed for the main analytic results 
of the article, which is given in the next section, for how to efficiently solve 
the Fermi-Bose model for a non-uniform molecular field. In particular we show 
how the system matrix for the dynamics of the Fermi-Bose model of a physical 
system in D spatial dimensions can be classified in terms of generalizations of 
Hankel matrices. 
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In cases where the shape of the molecular condensate p (x) posseses cer- 
tain geometrical symmetrices, such as spherical symmetry, a reduction of the 
number of atomic creation and annihilation operators can be implemented. Al- 
ternatively, a base with atomic operators that are directly defined e.g. in a 
spherical coordinate system can be used [35]. In this article, however, we treat 
the case of a general dense cubic lattice. In practice, the set of indices k (and 
k') lie on a finite grid in R D of the form k = (fci, . . . , ko) = x" 11 , where n £ IP 
and —K < rij < K for all j = 1, . . . , D. By abuse of notation, we will often 
suppress the factor 2jS and identify k with n, i.e. we will write e.g. ct nj i in place 
of ak,i wherever convenient. Whenever n or n' appears it will be implicitly 
understood that it stays within the above limitations. 



3.1 One-dimensional systems 



To set the scene, we first discuss a system in one spatial dimension. For a 
non-uniform system with D = 1, set B = 2K + 1 and identify the systems 
of annihilation operators {a nt i} -K<n<K and creation operators {a n %\—K<n<K 



with the 2B dimensional column vector 



-K,l ■ ■ ■ ttJf,l tt -A' 



K,2 



. Under 



this identification, Heisenberg equations (JS])-© can then be visualized in terms 
of a 2B x 2£?-system-matrix A composed of four blocks of B x B-matrices 



d 
~dt 



a-K,i 



aK,i 
t 

-K,2 



a 



L K,2 



A n A 12 
A 2 i A 22 



ax a 

t 

-K,2 



a 



l K,2 



(16) 



It then follows directly from © that the coupling matrices A12 and A21 become 
Hankel matrices. A Hankel matrix have the following structure A12 (m, n) = 
A\2 (m — 1, n + 1), i.e. the elements are identical when the sum of the row and 
column indices are constant |41) . 



3.2 Higher-dimensional systems 

For D > 1, the summation operator that takes the system of creation operators 
{ a k' 2}^' t° the system of annihilation operators {ak,i}k is defined via 

K K 

^12,n({aL, 2 }n') = <7 K H 9n 1 +n' 1 ......n D +n' D a i nli ,, n , DT (17) 

n' x — — K n' D — — K 

We call this a D-dimensional finite Hankel operator. These have been studied 
e.g. in |42| for D = 2. Obviously, there are multiple ways to represent the two 
systems {a^, 2 }k' and {ak,i}k as a 2B D — dimensional vector. However, with this 
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done, the Hankel structure of the corresponding coupling matrix A\2 is lost. In 
the next section we will construct such a concrete representation in which A12 
turns out to be a Z?-block-Hankel matrix. 



3.2.1 Ordering the lattice 

Recall that B = 2K + 1. The following function 

D 

f (m, n D ) = 1 + ( n i + K ) B °~ j > 
3=1 



(18) 



is a one-to-one mapping from the D-dimensional lattice n\, no, fij £ {— K, K} 
to {1, B D ). 
The inverse is 



«£> = mod (m — 1, B) — K 



f-i (m ) = J n d = mod (m - 1 - Ef= d +i ("i + *Q B fl ^,B fl+1 - d ) / B D ~ d - K. 



ni 



>D-1 



A" 



(19) 

Note that the equations (JTSJ) and (JTHJ) are really nothing else than a change of 
base from B to 10 for integer numbers. 

We will later also use the following property of the map above 



/- 1 (B D + 1 - m) = -r 1 (m) . 



(20) 



3.3 General construction of the system matrix 

We now identify the two lattices of operators {a n i} and {a^ 2 } with the corre- 



sponding column vector 



■ • ■ 0>B D S a \.2 ■ ■ ■ a B D 2 



, where we use the simpli- 



fied indices 1, B D instead of / 1 (1), / 1 (B D ). Then the D-dimensional 
system (H])-© can be represented in matrix form analogue to (fT6|) as 



d_ 

lit 



ai.i 



t 

1.2 



ft 

B D .2 



An A12 
A21 A22 



ai,i 



t 

1.2 



ft 

l B D ,2 



(21) 
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Clearly, An and A22 are diagonal matrices whereas A12 corresponds to the 
D-dimensional Hankel operator (fl7|) and A21 to the lower counterpart in ([9]). 
As the matrices, A12 and A21 are not Hankel matrices for D > 1, but rather 
exhibit a block-Hankel structure [32], we will refer to them as £>-block-Hankel 
matrices. In the coming two sections we discuss properties of the four block 
matrices An, Ai 2l A 2 i and A 2 2- 

We will use the following notation; T denotes transpose, * is complex conju- 
gation, and f represents both the two previous operations combined. Moreover, 
the operation of transposing in the skew-diagonal will be denoted SDT, i.e. 

A SDT (m R , m c ) = A(B D + 1 - m c , B D + 1 - m R ), (22) 

or equivalently 

A SDT = SA T S, (23) 

where S denotes the skew-diagonal identity, i.e. the matrix obtained by reversing 
the order of the columns of the identity matrix I. From Q23p and the property 
SS = I, it also follow that the skew-diagonal transpose of the product of two 
general matrices B\ and B2 fulfills 

(B 1 B 2 ) SDT = Bl DT B( DT , (24) 

which will be used later. Finally, the skew-diagonal transpose combined with 
complex conjugation will be denoted SDH . 

3.3.1 Structure of the Z?-block-Hankel matrices 

Inspection of ([H])-© an d ([21]) shows that A12 is given elementwise by 

A 12 (m R ,m c ) = qKg f -i {rn R )+f -i im cy (25) 

If we denote the coordinates of / _1 (m R ) by a sup-index R, those of / _1 (m ) 
with a C and the coordinates for the Fourier coefficients by , then the coor- 
dinates for the rows and columns of A12 fulfills 



(26) 



n D = n D + n D 



which considerably simplify any practical implementation of Ai2- Note also that 
when <7k is defined on the same k-lattice as the atomic operators, we necessarily 
have (jk = when \n R + > K, see figured] for an illustration. 
Similarly to ([23]) . we have 

A 21 (m R ,m c ) = ^-i (m « )+/ -i (mC) , (27) 

such that (q 2 = 1) 

A21 = qA* 12 - (28) 
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Due to this identity we are satisfied with discussing properties of A\ 2 in the 
remainder. First we note that 



A{ 2 = A 12 , (29) 

which is immediate by (|25p. 

To get additional structure, we impose extra conditions on the condensate 
wave-function "J = \J p (x) exp (id (x)). For many physical applications, "J = 
y/ p (x) can be chosen real. In this case, we note that the Fourier coefficients 
(fTU|) have the following symmetry g_k = <7 k . Therefore we get from (f2"0)) 

A^ H = A 12 , (30) 

and from (gS]) that 

A 21 =qAf 2 DT - (31) 

Furthermore, for the physically important case of a condensate wave-function 
that is real and even, i.e. with p(x) = p(— x), such as for example for a 
condensate in a harmonic trap, we have that <?k is real so A\ 2 = A\ 2 , which 
combined with (|3"0")) implies that 

Af 2 DT = A 12 , (32) 

and from (gS]) that 

A 21 = qA 12 . (33) 

Finally, note that in the real uniform case, p(x) = po, = 0, we have non- 
zero entries only at positions where (n^, n^J = — (nf , n^) , i.e. Ai 2 = 
qA 2 i = qngoS, as in section T2.2I 



3.3.2 Diagonal matrices for the kinetic energy 

According to the Heisenberg equations (j8])-(j9]) and the representation (|2"Tj). the 
matrices An and A 22 become diagonal (kinetic energy in momentum space), 
where the m:th diagonal element of An is given by 

A 11 (m,m) = -iA k = -iUl + ^-Y m = f(k), (34) 

and similarly 

A 22 (m, m) = iA k = A* n (m, m). (35) 

Finally, we note that 

Afr = An, (36) 



which follows directly from (|20j) . since |k| = | — k 
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4 On the block-structure of the propagator 



Much of the structure of the system-matrix A is preserved in matrix functions 
defined on A, which can be used to reduce the computational complexity when 
evaluating exp (At) for the system's evolution in time. Naturally, the more con- 
ditions we impose on the condensate wave- function in (fTOj) . the more struc- 
ture is preserved. We present the corresponding identities in order of increasing 
symmetry on vff, starting with a general complex condensate wave-function. 
Throughout we will let q = ±1 and prove the results for the cases of fermionic 
(q = —1) and bosonic (q — 1) atoms simultaneously. 

Given an arbitrary square (even sized) 2n x 2n-matrix B we decompose it 
into its four blocks 

B\\ B12 
B21 B22 



B 



Below we list a number of useful matrix identities which all can be verified by 
direct computation 







B21 B22 



I 
ql 



B22 

qB\2 



qB 2 i 

B u 



(37) 



B21 



B\2 

B22 



Bl 
Bl Bl 



(38) 



Bu 
B21 



B\2 

B22 



SDT 



B SDT 
D 22 
B SDT 
n 21 



B SDT 
n 12 
B SDT 
D U 



(39) 



With these three identities at hand, the following lemmas are all immediate. 
For example the first one is a direct application of (|3"T)) alone. 

4.1 Lemma 

The matrix identities 



are equivalent to 



Bu — B22, B12 — qB^i, 



B 



ql 

1 



B* 



I 
ql 



(40) 
(41) 



The second lemma follows by combining (|57|) with . 
4.2 Lemma 

The matrix identities 



r> tjSDT 72 r>SDT r> 

&U — &u ; ^22 — £>22 J "12 



r>SDT 
QB 21 , 



(42) 
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are equivalent to 



D 



D 



SDT 



I 
ql 



(43) 



Finally combining (|3~T|) with (|3"5)) we obtain the third lemma. 



4.3 Lemma 

The matrix identities 



Bn — B\ 2 , B12 — b\ 2 , B 2 \ — B\ x , 



are equivalent to 



B 



I 

1 



I 

1 



(44) 
(45) 



4.4 Proposition 

Let B be the system-matrix A constructed in section \ 3.3\ from the condensate 
wave-function ^ . Then the identities in j^fl[ ) of Lemma \4--1\ are always satisfied. 
Moreover, if ^ is real then the identities in J^[ ) of Lemma \4-2\ hold and if ^ is 
real and even, the identities in \44\) °f Lemma \4-3\ hold. 

Proof. The first claim follows by combining (|28|) and (|3"5"T) . and the second 
follows by combining (|3"Tj) and ([3T>|) . When "F is also even the elements of A 12 
are real, and hence (|29[) can be written A\ 2 = A\2 (analogously A\ x = A21 ). 
Since clearly An — A 22 , the last claim is established as well. 

□ 

We need some preparation for the main theorem and corollaries, which basi- 
cally says that the above identities are preserved when forming exp (At). First 
we note that 






I " 




' 


ql - 




' I 




I " 


. I 1 







I 







q 2 l _ 




I 



(46) 



We now recall that an analytic function <j>(z) = J2T=o Ckzk defined on all of 
C is called an entire function, and that for each convergence radius r > one 
can find a constant C r such that |cfe| < [43]. This allows us to define the 
function 4>(B) for any matrix B (with the corresponding matrix norm ||-B||), 
since 

N 2 N 2 N 2 11 

II £ rfll < £ kii|B fc || < £ Cr^-S-, 

k=N 1 k=N! fc=iVi 

which (picking any r > ||-B||) shows us that (X^fclo c kB k )^ =0 is a Cauchy se- 
quence in the set of matrices {B} with the operator norm ||-B||. The limit is 
thus a well defined matrix since this is a Banach space. We recall that ||f?|| is 
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defined as \\B\\ = sup{||B(x)|| : ||x|| = 1} and in particular that for any given 
index (to, n) we have 

\B(m,n)\<\\B\\. (47) 
We are now ready for the main result. 



4.5 Theorem 



Let 4>(z) — X)fco c k zk be an entire function and let A be a matrix that satisfies 
either of the identities in \4'/fy or \44ty - Then <fi(A) satisfies the same 

identities. 

Proof. Let us suppose that A satisfies the identities of (|4"!2)) . By Lemma [ 
we then have (|43p . which combined with (|24p and (|4l)|). yields that 



A k 



ql 

1 



(A 



k\SDT 



I 

ql 



for any k. Thus the matrix identities in (|42[) are satisfied for the corresponding 
blocks of M = 4>{A) when is a monomial. Since (|4"2j) are also preserved upon 
taking linear combinations of matrices that all satisfy (|4"2)) , we conclude that 
(|4^|) holds whenever <j> is a polynomial. Finally, since in the general case we 
have 

N 



cj>{A) = lim Vc^, 



k=0 

in the operator norm, and since the identities (|42p are preserved upon taking 
limits with respect to this norm [see (|47|) ]. the general result follows. The proofs 
related to (|4"0|) and ((U]) are analogue. 
□ 

Remark: Connecting to Proposition 14.41 we note that when B arise from 
a general complex condensate wave-function, one still has more structure than 
what is expressed in (|4"U|) . For example 

— Bn = B T 2 , B 12 = Bj 2 , B 2 \ = B 21 , 

which is easily seen to be equivalent to 



B = 



-I 

1 



B 1 



-I 

1 



' -I " 




— I " 




-I 


I 




I 




-I 



However, since 



these properties are not preserved in 4>{B), except for in the special case when 
4> is an odd function. 

We now sum up our conclusions for the physically most interesting cases. 
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4.6 Corollary 

Given a real condensate wave-function ^ and f G 1, set 

M = exp(^)=E^- 



Then AI has the structure 



M 



M n 
M* 2 



3=0 



qM 12 
M u 



where in addition we have for the two blocks 



(48) 



(49) 



7, rSDT 



M n , M? 2 DH 



Mi 



(50) 



Proof. We can see that M has the above structure if and only if it satisfies 
the identities in (j40|) and (|42j) of Lemmas 14.11 and 14.21 Moreover, A satisfies 
these identities by Proposition 14.41 The desired conclusion thus follows by The- 
orem 14.51 

□ 

It is clear from Corollary 14.61 that we only need to calculate half of the 
matrices Mn and M 12 in order to fully determine M, which generally reduces 
the computational cost from (2n) = An 2 to 2 (n 2 /2) = n 2 elements in this case. 



4.7 Corollary 

Suppose that ^ is a real and even condensate wave- function. Then, in addition 
to the identities in Corollary \4-.6[ we have 

Mfj = Mn, M 12 = Mi 2 . (51) 

Proof. By Proposition 14.41 A satisfies the identities in (|4"4"|) of Lemma 14.31 
and hence so does M by Theorem 14.51 It is easy to see that these identities 
combined with the structure in (|49|) and (|50"|) proven in Corollary 14.61 implies 
that the above identities are satisfied for M . 

□ 

From Corollary 14.71 follows that we only need to calculate a quarter of the 
matrices Mn and Mi 2 in order to fully determine M, which further reduces the 
computational cost to 2 (n 2 /4) = n 2 /2 elements in this case. 



5 Obtaining physical observables 

In this section we show how to use the results of the previous section in obtaining 
physical observables for the atoms. We start from the following general block 
form of the solution M = cxp (At) to Heisenberg equations ©-([S]) in matrix 
form (IT]) 
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oi,i (*) 



a B D A (t) 
* :(*) 



'1,2 



l B D ,2 



(*) 



Mn 0M12 
M 2 i M 22 



«M (0) 



a B n A (0) 
*t (0) 



"1,2 



it 

i B r> .2 



(0) 



(52) 



It is obvious that the results of the previous section will simplify the practical 
calculations of M (t) and hence the physical observables. However, first we 
show in the next section how to generally obtain first-order moments for pairs 
of atomic operators directly from (|52[) . 



5.1 First-order atomic moments 

We now denote by My^ the m-row- vector of the block matrix My, where 
m = f (n) is mapped to = k = ki, ko according to section l3.2.1l 

As a first example, for an annihilation operator of the a = 1 spin-state in 
row m in the left hand side of (|52p we have 



a m ,i (*) = o k ,i (t) = Mu,kU + qMi 2 ^v = u T k + qv T M^ 2 k , 



(53) 



where we for notational and computational convenience introduce the two op- 
erators 



(0) 

3 B D 1 (0) 



v = [ 01,2(0), . . . ,050,2(0) 



4,2 (0) 

a BD (0) 



which are naturally constructed from (|52[) . With u and u we have introduced 
in effect a calculus, where only terms in the expectation values containing the 
matrix (uv)^ = I or (i^ v T ^ = I give a contribution, while all other terms 

are zero. This is due to the (anti-) commutator relations applied to an initial 
vacuum state for the atoms. See the examples leading to ([55]) and ([57)) below for 
further details. With these rules, different vacuum expectation values of atomic 
operator pairs can conveniently be written down in terms of standard matrix 
products between complex row- and column- vectors, where both the vectors are 
defined from certain rows in one of the four block matrices . 

We now give the corresponding expression to (|53[) for an annihilation oper- 
ator of the er = 2 spin-state 



Ok,2 (*) = (o 



"k,2 



(M 2 i,kM + M 2 2, k 6) 1 



fitMt 



(54) 
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Note that analogously to ([55)1 , we can also alternatively write (fMf on the form 

Qk. 2 (t) = M£i + ^22 i an< ^ f° r a specific moment one chose the option 
that allow the operators to meet in the middle of the operator pair. 

From (|5U)) and ([5^)1 we can then demonstrate in detail how to calculate the 
anomalous moments, 



mk.k- = (ak.iSk'.a) = |(Mn,k« + gMi2,k«) f^-^li.k' + f,t -^22,k' 
Mn ik uu t M 2 t lk , + gM^.k^+M^k, 



+Mu, k m) t M 2 t 2 k , + gMi2,kU* + A4 2jk , 
= Mu ik (rf) M| lik , + gMi2, k (w^) A4 ljk , 

+ M u , k (tifit) M 2 f 2 k , + gA/ 12 , k (««t) M 2 f 2)k , = M llik Mj 1|k „ (55) 

which is a time-dependent complex number as expected. 
In a similar way, using the Hermitian conjugate of (|53|) . 



4 i = (^n.kM + gA/i2,kw) f = u^-^i k + qvU4 2 k = A/*! k w tT + gM£, k ?) tT , 

(56) 

we have for the normal moments of the a=l spin-state 



\qM* l2 ^ v T qM? 2M ,) = q 2 M* 2M M? 2M , = M^M^,. (57) 
From ([Ml) it follows that the corresponding result to (|5"7| for the er=2 spin-state 



is 



"k,k',2 = (a k 2 a k %2^ = ^(Af2i,kW + M22,k«) f'" t M 2 t 1;k , + &M% 2 k , 

[M 3 i,k«« + -a4 lik ^ = M 21M M\ 1M ,. (58) 
We finally confirm by direct calculations from and (|56l) that 



™k,k' = ((ak.iak' : 2) f ) = (S k , i2 a k ,i) = M 21M' M Ii,\l = ( M n,k^ 2 Vk') > 
as expected from (|55l . 
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Any other first-order moment that can be formed by two atomic operators, 
such as for example 1 ak',2^i can be shown to be zero. 
We now note that with k' = k in (IST1) we have 



™k,l = ™k,k,l = \ a k,l fl k,l ) = ^12,k^l2,k 



which is a real number as expected from the physical interpretation of the 
occupation of spin-1 atoms in state k. 

Motivated by ([57)1 and (|58p we see that the fact that the two spin-states 
are analogue in the Fermi-Bose model ([T]), with initial atomic vacuum in both 
spin-states, physically supports the identity M21 = Mf 2 given in Corollary 14.61 
Hence, applying Corollary 14.61 it follows from (|5"7[) and ([55)1 that we can write 
for arbitrary spin 



(59) 



= \ a M fl k',"/ = M 12,k M 12,k' 

Similarly the anomalous moments from (|55[) then take the form 



m KW = M u k M 12 k ,, 



(60) 



in terms of only the blocks M\\ and M12. 



5.2 Deduction of the uniform case 



In this section we will deduce the solutions (fT3|) and (|14|) direct from the expo- 
nentialmatrix in order to check the formalism presented. We start by rewriting 
according to 



M 



exp (At) = £ 



£ 



(2j) 



£ 



(2J + 1)! 



(61) 



3=0 J j=o v ■" j=0 
Then starting from the system matrix for a uniform molecular field with go 
K 9o = Xy/PO 



A 



diag (^-iA k 
5o§ 



diag (iAkj 

we first note, using § 2 = I, that the even powers of A are diagonal, hence we 
have 



E 



(At) 



2/ 



^ (2i)! 
i=o v ■" 



diag cos y A k - qgfi t 



Secondly, we see by writing the odd powers as AA 2 i that 



diag cos v/ A k - qgfi t 



(62) 
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Finally, we can from (fSTj) . (j62|) and (|63|) write up M on the form of (|52|) 




From (64]) we can then via (J59J) and (JeUJ) retrieve (13]) and dHJ) of section l2~2l 



However, it is important to stress that the (skew-) diagonal blockmatrices of 
M lose this structure for a non-uniform condensate wave-function, and in the 
general case it follows by Cauchy-Schwarz that the equality (fl3|) changes into 
the following inequality 

|wk| 2 < n^ a (1 + qnk, a ) ■ (65) 

The above inequality have been investigated numerically for the case of disso- 
ciation into bosonic atoms in |24| and it implies limitations on the strength of 
various correlations for non-uniform systems |23[ 124] . 



5.3 Higher-order atomic moments 

Higher-order moments, such as in the simplest case, the combination of two pairs 
of operators, are factorized according to Wick's theorem [33] which is implicit 
from the decorrelation assumption in use within the undepleted molecular field 
approximation here. As an example we calculate Glauber's correlation function 
for two atoms in the same spin-state |23| 



g W ( k , k', t) ee \ k ' CTk ' g / = 1 + gl " k *' CTl ■ (66) 

"k.crnk'.cr «k,k,<r^k',k',(T 

For a numerical implementation of (|66[) we have, according to ([J 
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g™ (k,k',t) = 1 + ? - T? - T . (67) 

We here also give Glauber's correlation function for two atoms in opposite spin- 
states 

foi ("k l"k' 2"k',2«k.l) \ m , ,,| 2 

fl ? (k,k',i) ee \ k ' lk ' 2 / = 1 + |mk ' k '' ■ (68) 

rik.ink',2 «k,k,i«k',k'.2 

For a numerical implementation of (|68|) we have, according to ([59| and f|60f> . 





Mu )k A^ ik , 


2 


(^r 2) k^l T 2,k; 


(M* 2k ,Mf 2k/ ) 



5 J 2 2 ) (k,k',t) = l + 7 r . (69) 



5.4 Calculations of exponential matrices in practise 

Up to this point, we have shown how to reduce the computational needs for 
obtaining physical observables to the calculation of only a fraction of the full 
exponential matrix AI = exp (At) . We now discuss how we perform the neces- 
sary numerical calculations for the physical observables in practice, while the 
results are presented in the next section. 

The matrix exp (At) can be calculated numerically in many different ways, 
see for example [U] for a review of methods. For example, the two most obvious 
ones are: i) use the definition in (|4"5|) and calculate all powers until some cut- 
off in k, (At) k jk\ m 0, or = if A is nilpotent; ii) diagonalize A and use all 
its eigenvalues (Xj) and eigenvectors (sj) to change to the S'-basis, such that 
AI = SDS' 1 , where Djj = exp(Ajt). Both the natural methods i) and ii) 
in practice needs many matrix-matrix operations, and have computationally 
expensive performance for large matrices. 

As explained in section 15. 1[ all observables are obtained as matrix products 
between different row- and column-vectors defined from rows of the blocks of 
AI. It is obvious that row R of M is obtained by multiplying M T with the unit 
vector e# = [0, 1^, 0] T from the right. For this task, powerful matrix- 
free algorithms exists for general matrices [45], i.e. that calculate the results 
of an exponential matrix acting on an arbitrary vector, without the need to 
perform any matrix-matrix operations. 

Using the results of the Corollaries in section|4j the fact that An is diagonal 
and hence can be represented as a vector, and the powerful property of A12 
being a Z?-block-Hankel matrix, we use our own modified version of the open 
access software Expokit for sparse matrices [45]. As will be reported elsewhere, 
we have optimized the algorithm from |45| for the implementation of D-block- 
Hankel matrices. This software optimization, combined with the use of the 
results in section 21 give us the crucial advantages necessary for implementating 
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large Fermi-Bose systems, compared to any brute force calculation of the full 
exponential matrix. While the complex matrix M, or any of its blocks, is 
generally not sparse, the matrix A12 is sparse especially after truncation of the 
smallest Fourier coefficients, see figure Q] of the next session for an example. 
Clearly any such actual truncation have to be evaluated with convergence tests. 

Overall, our optimized exponentiation procedure in effect reduces the orig- 
inal problem of calculating the complex (non-sparse) 2n x 2n-matrix M = 
exp (At) to the simpler problem of performing matrix-vector operations with 
the D-block-Hankel matrix Aj 2 = A12, which in addition obey further sym- 
metries, see e.g. (f2U|) and ([50)) . Furthermore, A12 is real for the physically 
important case of a condensate wave-function that is even, which reduces the 
information to store in A12 by an additional factor of two in this case. 

Finally, let us stress the important consequence of section [5Tl that the cal- 
culation of each row in the blocks My is independent, such that the work can 
conviniently be distributed across several computers in parallel to reduce com- 
putation time if needed. This is in principle a very crucial advantage of the 
presented formalism, when applied to large systems. 

6 Numerical illustrations 

We here present specific numerical results for atomic correlation functions for a 
non-isotropic three-dimensional system. We use a 61 x 61 x 61 grid in momentum 
space, which corresponds to a system of linear differential operator equations 
with a 2n x 2n-system matrix where n = 61 3 = 226981, i.e. in general (2n) 2 ~ 
2.1 • 10 11 matrix elements. The theorem presented in section 2] in combination 
with the algorithm improvements briefly discussed in section !?^! allow systems of 
this size to be solved on one standard PC in the order of ~ 10 hour. As pointed 
out in section HOI substantially larger grids can be attacked with calculations 
on several computers in parallell. 

6.1 Physical parameters 

In order to use a realistic set of parameters for our numerical example, we 
chose an harmonic trap with frequences such that the size of the molecular 
field along the x-direction Rtf.x = 8 /im is two times the size along the z- 
direction Rtf,z = 4 /mi, while the size along the y-direction RrF,y = 6 /<m is 
set to an intermediate value. With a central molecular peak density of po = 
10 20 m -3 , this corresponds to No ~ PoRtf,xRtf,i/Rtf,z — 3.2- 10 4 molecules. 
Choosing 40 dimers [13] we have an atomic mass of m at = 6.642 • 10 -26 kg. 
The molecule-atom dissociation parameter is x — 10 _7 m 3 / 2 /s here |2H]. We set 
the dissociation detuning to ft = — 4- 10 3 s~ 1 which is large enough to ensure that 
the dissociation energy is larger than thermal excitations at nK temperatures 
(i.e. 27i|f2| 3> ksT). With a characteristic time of to = 1ms this corresponds 
to a dimensionless detuning of 5 = to^l = —4 }28j. The momentum lattice in 
use have a spacing dk = dk x = dk y = dk z ~ 1.1 ■ 10 5 m _1 which is smaller than 
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Figure 1: Illustration of the D- block- Hankel structure for the actual A12 block- 
matrix used in the numerical example with D = 3 and K = 30 reported in this 
section. Figure (a) illustrates the entire 226981 x 226981 A12 matrix, while (b), 
(c) and (d) shows zoomed in regions corresponding to the circle of the previous 
subfigure. Note that in (b) each dot represents a 2-block-Hankel matrix, while in 
(c) and (d) each dot represent a single non-zero matrix element. Dashed (dotted) 
lines in (a) shows the skew-diagonal (diagonal) which corresponds to the skew- 
diagonal transpose (SDT) (respectivelly SDH for complex matrix elements) and 
transpose symmetries discussed in section [3.3. II However, let us stress that the 
symmetry of a D-block-Hankel matrix is considerable higher, since the block- 
matrices in (a) and (b) [elements in (d)] are repeated downwards parallell to the 
respective skew-diagonals. Grey triangular regions in (a) and (d) shows where 
the conditions |nj* + > K for j = 3 respectivelly j = 1 are fulfilled. As 
discussed in section 13.3.11 these regions does not contain any non-zero matrix 
elements due to the restricted Fourier coefficients lattice size. The fact that the 
regions of non-zero elements do not reach out to the grey fields are due to the 
truncation implemented, see text for details. 
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the smallest width of the molecular momentum distribution ~ 2/Rtf,x 123) , 
and have been confirmed numerically to resolve the dynamics of the relevant 
structures in the atomic momentum distribution. The corresponding resonance 
momenta is then fco = |ko| = ^/2m a t [fij /H ~ 20<ifc. 

6.1.1 Fermi's Golden Rule estimate of the atom numbers 

For small times we can, for the purpose of validating the physical parameters, 
estimate the number of atoms by the following linear expression in time 

Nj (t) c iV At, A = -A- (^f'\ 2 ^T\. (70) 

From the above formula it is clear that the number of atoms increase with |f2| 
for a three-dimensional system. This have earlier been studied explicitly for 
a uniform system, see Fig. 1 of |28| . Hence for cases with large detuning the 
validity of the undepleted field approximation is limited to short dissociation 
times t/to <C 1. For the parameters of section we have A ~ 2s _1 which 
results in Nj ~ 10 2 atoms at t = to, and hence a conversion ratio of less than 
1%, ensuring the validity of the results from the undepleted field approximation 
[23l [25l [26] . The presented estimate of atom numbers from the Fermi's Golden 
Rule (|70| was later confirmed by the numerical calculations for the parameter 
values in use here. 

6.2 Structures in the system matrix 

We here explicitly illustrate the Z?-block-Hankel matrix A12 that is used in the 
numerical calculations of a physical system for D = 3 here. Hence, it is evident 
from the figures [1] (a), (b) and (c) that we can zoom in D = 3 times on A\2 and 
reveal a repeating pattern. Clearly after the last zoom in [figure [1] (d)], we are 
left with a structure of a usual Hankel matrix. 

We have used a truncation of the molecular BEC source such that Fourier 
coefficients with a modulus less than 2% of the leading coefficient is neglected. 
This procedure have been evaluated by reconstruction of the BEC by the inverse 
Fourier transform, ft was also found that the 2% level of truncation resulted in 
correlation functions (see figure ^ that could not be distingushable by the eye 
from the correlation functions obtained with a 4% level of truncation. 

6.3 Numerical comparison with analytic asymptotes 

For D = 3 the Thomas-Fermi (TF) density profile of the molcular BEC is 
given by p {x) = p {\ - x 2 jR\ Fx - y 2 /R\ F , y - z 2 /R 2 ?Fz ) for x 2 /R 2 :Fx + 
y 2 J 1 R\ F y + z 2 /R 2 ?F z < 1 [and po( x ) = otherwise], which is underlying an 
analytic derivation of the asymptotes. Here Rtfj is the Thomas-Fermi radius 
along the spatial direction j — x,y,z. We are here interested in collinear (CL) 
density correlations between two momentum components at k and k' , for which 
the displacement Ak = k k' is along one of the Cartesian coordinates, kj. The 
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2 t/t 2 tft 

k /k„ k /k„ 

x z 

Figure 2: Fermionic collinear atom-atom correlation functions in momentum 
space gjy (k, k , t), at times t/to = 0.1, 0.2, 1 (to = 1ms), calculated along dif- 
ferent directions in 3D. In (a) we show numerical results (solid thin curves) from 
(|67j) along the direction e^, i.e. with k = fc^e^ and k' = koe x (fco — 2.2-10 6 m~ 1 ), 
while in (b) we show the corresponding result along the direction e z . The ana- 
lytic short-time asymptotes of (J7TJ) are represented by the fat curves plotted at 
t/to = 0.1 only. In fact the short-time asymptotes for the collinear correlations 
are in qualitative agreement with the numerical results even up to t/to ~ 1- 
However, a zoom in reveals quantitative deviations seen as a narrowing of the 
width of the correlation signal with time, this is also in agreement with detailed 
ID results reported in Fig. 8 of [53]. In general the fermionic collinear correla- 
tions are here showing a Pauli-blocking dip at k XyZ = fco, while the characteristic 
width of the correlation signal have been confirmed to be inverselly proportional 
to the size of the molecular BEC source along the corresponding direction, i.e. 
- 2.l6R^p x ~ 2.7 • lOW 1 in (a) and - 2.16-R^ = 2 • 2.16i? T j ?:c in (b). In 
addition we compared with the corresponding results for bosonic atoms showing 
a so called Hanbury-Brown and Twiss peak at k x , z = (dashed curves), shown 
only for the largest time here. 
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detailed derivation of short-time asymptotes for the correlation functions in this 
case was reported in [23] . The CL correlations following from this derivation is 

9 ^,^l + ^ lJ '^ k '-^f ] \ (71) 

where J„ denotes Bessel functions of the first kind. The qualitative behavior of 
the CL correlation functions are similar as in lower dimensions |34( 146j , whereas 
the quantitative differences enter e.g. through the width and the peak values. 
The widths of (|7Tj) is w\ ' ~ 2.16/RrF.j and the peak value of (|71|) is static 
to leading order |23) . 

In figure [2] we show results for the evaluation of the analytic short-time 
asymptote (fTTjl . strictly valid in the t/to <C 1 limit, against numerical results 
for times t/to < 1, where t = to roughly corresponds to the first maximum in 
time of the oscillating fermionic atom numbers N (t) = ^ k rik.er {t). 



7 Summary 

We have described how to effectivelly calculate the dynamics of linear Heisen- 
berg operator equations for the Fermi-Bose model applied to the problem of 
molecular dissociation. We note that a similar framework have been used to 
obtain numerical results for a non-isotropic 2D system on a 61 x 61 grid in 
|34) . We have here generalized the approch to D spatial dimensions with the 
use of D-block-Hankel matrices. In particular we have explicitly explored a 
non-isotropic 3D system on a 61 x 61 x 61 grid numerically on a standard PC. 
Such a grid can resolve relevant atom dynamics in momentum space for realis- 
tic parameters [9], and naturally extends earlier studies of non- uniform ID and 
2D systems [341 146j , and is more realistic than previous treatments of uniform 
3D systems [HI [2E] ■ We finally stress that the results presented can be used to 
handle a complex bosonic mean-field of any geometry in any spatial dimension. 

Acknowledgments 

We thank Karen Kheruntsyan and Roger Sidje for valuable discussions at an 
early stage, and Johnny Kvistholm for artistic assistance with figure 1. 

References 

[1] Friedberg R and Lee T D 1989 Phys. Rev. B 40, 6745. 

[2] Kheruntsyan K V and Drummond P D 2000 Phys. Rev. A 61, 063816. 

[3] Holland M, Kokkelmans S J J M F, Chiofalo M L and Walser R 2001 Phys. 
Rev. Lett. 87, 120406. 



26 



[4] Timmermans E, Furuya K, Milonni P W and Kerman A K 2001 Phys. Lett. 
A 285, 228. 

[5] Ohashi Y and Grin A 2002 Phys. Rev. Lett. 89, 130402. 

[6] Strohmaier N, Greif D, Jordens R, Tarruell L, Moritz H, Esslinger T, Sen- 
sarma R, Pekker D, Altman E and Dernier E 2010 Phys. Rev. Lett. 104, 
080401. 

[7] Jordens R, Strohmaier N, Gunter K, Moritz H and Esslinger T 2008 Nature 
(London) 455, 204. 

[8] Jack M W and Pu H 2005 Phys. Rev. A 72, 063625. 

[9] Kheruntsyan K V 2006 Phys. Rev. Lett. 96, 110401. 

[10] Mukaiyama T, Abo-Shaeer J R, Xu K, Chin J K and Ketterle W 2004 
Phys. Rev. Lett. 92, 180402. 

[11] Diirr S, Volz T and Rempe G 2004 Phys. Rev. A 70, 031601 (R). 

[12] Thompson S T, Hodby E, Wieman C E 2005 Phys. Rev. Lett. 94, 020401. 

[13] Greiner M, Regal C A, Stewart J T and Jin D S 2005 Phys. Rev. Lett. 94, 
110401. 

[14] Lee T D 1954 Phys. Rev. 95, 1329. 

[15] Kallen G and Pauli W 1955 Kgl. Danske Vidensk. Selsk. Mat.-Fys. Medd. 
30 No. 7. 

[16] Heisenberg W 1957 Nuclear Physics 4, 532. 

[17] Poulsen U V and M0lmer K 2001 Phys. Rev. A 63, 023604. 

[18] Drummond P D and Kheruntsyan K V 2002 Phys. Rev. A 66, 031602(R). 

[19] Cazalilla M A and Marston J B 2002 Phys. Rev. Lett. 88, 256403. 

[20] Schollwock U 2011 Annals of Physics 326, 96. 

[21] Deuar P, Chwedenczuk J, Trippenbach M and Zih P 2011 Phys. Rev. A 83, 
063625. 

[22] Krachmalnicoff V, Jaskula J-C, Bonneau M, Leung V, Partridge G B, Bo- 
iron D, Westbrook C I, Deuar P, Zih P, Trippenbach M and Kheruntsyan 
K V 2010 Phys. Rev. Lett. 104, 150402. 

[23] Ogren M and Kheruntsyan K V 2010 Phys. Rev. A 82, 013641. 

[24] Savage C M, Schwenn P E and Kheruntsyan K V 2006 Phys. Rev. A 74, 
033620. 



27 



[25] Midgley S L W, Wiister S, Olsen M K, Davis M J, and Kheruntsyan K V 
2009 Phys. Rev. A 79, 053632. 

[26] Ogren M, Kheruntsyan K V and Corney J F 2010 Europhys. Lett., 92 
36003. 

[27] Ogren M. Kheruntsyan K V and Corney J F 2011 Comput. Phys. Commun. 
182, 1999. 

[28] Davis M J, Thwaite S J, Olsen M K and Kheruntsyan K V 2008 Phys. Rev. 
A 77, 023617. 

[29] Corney J F and Drummond P D 2004 Phys. Rev. Lett. 93, 260401. 

[30] Corney J F and Drummond P D 2006 J. Phys. A: Math. Gen. 39, 269. 

[31] Rahav S and Mukamel S 2009 Phys. Rev. B 79, 165103. 

[32] Rosales-Zarate L E C and Drummond P D 2011 Phys. Rev. A 84, 042114. 

[33] Fetter A L and Walecka J D 2003 Quantum Theory of Many-Particle Sys- 
tems, Dover. 

[34] Ogren M, Savage C M, and Kheruntsyan K V 2009 Phys. Rev. A 79, 043624. 

[35] Deuar P, Zifi P, Chwedehczuk J, Trippenbach M 2011 Eur. Phys. J. D 65, 
19. 

[36] Poulsen U V and M0lmer K 2007 Phys. Rev. A 76, 013614. 

[37] Castin Y and Dum R 1996 Phys. Rev. Lett. 77, 5315. 

[38) Zih P, Chwedehczuk J, Trippenbach M 2006 Phys. Rev. A 73, 033602. 

[39] Ogren M and Kheruntsyan K V 2009 Phys. Rev. A 79, 021606(R). 

[40] Ogren M and Kheruntsyan K V 2007 unpublished. 

[41] Peller V 2003 Hankel operators and their applications, Springer- 
Verlag N. Y. 

[42] Andersson F, Carlsson M and de Hoop M V 2010 Appl. Comput. Harmon. 
Anal. 29, 156. 

[43] Lang S 1999 Complex Analysis, 4.th ed. Springer- Verlag N. Y., Grad. Texts 
in Math. 103. 

[44] Moler C and Van Loan C 2003 SI AM Rev. 45 3. 

[45] Sidje R B 1998 ACM Trans. Math. Softw., 24(1):130-156. 

[46] Ogren M and Kheruntsyan K V 2008 Phys. Rev. A 78, 011602(R). 



28 



